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Abstract 

The purpose of the project was to investigate methods to accurately verify that designed 
materials meet thermal specifications. The project involved heat transfer calculations and 
optimization studies, and no laboratory experiments were performed. One part of the 
research involved study of materials in which conduction heat transfer predominates. 
Results include techniques to choose among several experimental designs, and protocols 
for determining the optimum experimental conditions for determination of thermal 
properties. Metal foam materials were also studied in which both conduction and 
radiation heat transfer are present. Results of this work include procedures to optimize 
the design of experiments to accurately measure both conductive and radiative thermal 
properties. Detailed results in the form of three journal papers have been appended to 
this report. 


Introduction 

Research on functionally graded materials has been underway for several years. 
Functionally graded materials include composites with epoxy and metal matrices; metal 
foams; and any built-up material or structure with properties designed to vary across the 
structure. In the future, as functionally graded materials are incorporated into aerospace 
structures, part of the procurement process will involve certification that the materials 
meet the specifications. However, there has been little research on accurate thermal 
characterization of functionally graded materials. The present project is a first step 
towards closing this gap in the procurement cycle by developing accurate methods to 
measure the thermal properties of functionally graded (FG) materials. 

Two classes of materials were studied in this project. First, materials with spatial 
variation of thermal conductivity and specific heat were studied. If the maximum 
(simulated) temperature is kept low, the primary mechanism for heat transfer is 
conduction (molecule-to-molecule diffusion). Simulated experiments were studied that 
involved transient heating, and also steady-periodic heating. Second, metal foams were 
studied for which both conduction and radiation are important in the transfer of heat. 
Metal-foam materials are being studied as possible components of aerospace thermal 
protection systems (TPS). 

In this report, a brief account of the results of each of three phases of the study is 
given. Detailed results are given in the appendix in the form of research papers. 


Conduction-dominated material, transient heating 

The first phase of the research involved transient-heating experiments for the 
measurement of thermal properties in FG materials with spatially-varying thermal 
properties. The maximum (simulated) temperature rise was limited to about 40 C so that 
conduction would be the dominant mechanism for heat transfer. 

One-dimensional transient experiments were studied for a material with linearly- 

varying thermal properties described by parameters C , k , and spatial slope e. Variation 
of properties with temperature was not treated. Several transient experimental designs 
were considered, including single or multiple heating events and including surface and/or 
interior temperature sensors. An optimality criterion D, based on sensitivity coefficients 
[1] was used to find the best operating conditions for each experiment, and to find the 
best experimental design among those studied. 

The results of the optimality study show that it is more difficult to obtain accurate 
parameters when spatial slope e is small compared to when e is large; that is, it is easier 
to "see" large spatial variations in properties. The best experimental design involves 
analysis of combined data from two separate heating events, one with heating on one side 
of the body and one with heating on the other side, each time with the sensor located on 
the heated side and the unheated side maintained at a fixed temperature. This conclusion 
is also supported by a series of simulated experiments, carried out with regression 
analysis of error-containing temperature values. The optimal operating conditions for the 
experiment are somewhat dependent on the property slope e. Detailed results are given 
in Appendix A. 


Conduction-dominated material, steady-periodic heating 

This phase of the research investigated steady-periodic methods for thermal 
characterization of functionally graded materials. Specifically, the work was carried out 
to explore, via numerical simulation, the use of photothermal methods. A new heat 
transfer model was developed for the response of layered material to periodic heating, 
based on the method of Green's functions, which is numerically better behaved than 
previous work. The method was applied to a Si02 layer on silicon and compared to 
literature values [2] to validate the method. The optimality criterion D, introduced earlier, 
indicates which frequency range of experimental data should provide the best possible 
estimates of the layer conductivity and contact conductance. 

The new steady-periodic model was then applied to a two-component 
functionally-graded material with a power-law distribution of thermal properties. This 
type of material was studied previously for its thermal-stress behavior [3]. The results 
indicate that largest temperature response is found by heating the sample on the low-& 
side, and the phase of the temperature is particularly important for estimation of thermal 
properties. Component conductivities ki and k 2 have similar-shaped sensitivity 
coefficients, and consequently both cannot be estimated simultaneously from 
experimental data. The most important parameter is power-law exponent p which 
describes the spatial distribution of thermal properties in the functionally-graded material. 
Values for optimality criterion D indicate that values for p may be found simultaneously 
with one of the conductivities, but not both. Dimensionless frequencies less than unity 
are important for measurement of spatial-distribution parameter p. The magnitude of the 
optimality criterion D also suggests that it will be easier to estimate parameters for p - 1 
(near-linear spatial variation) compared to other values of p. More detailed results are 
given in Appendix B. 


Metal Foams with conduction and radiation heat transfer 

In this phase of the research the methods of optimal experiment design are applied 
to a high-porosity nickel foam material, with thermal property values taken from Sullins 
and Daryabeigi [4], Contained within the model are five parameters considered to be 
intrinsic properties of the material that must be determined empirically. Two of these 
intrinsic properties are related to heat transfer by conduction and three are related to the 
radiative heat transfer. The intrinsic properties are the following: F, the efficiency for 
solid (metal) conductivity; a, coupling coefficient for solid-gas conductivity; eo and e u 
parameters defining the temperature-dependent specific extinction coefficient as 
e = eo + ej T; and, co , the albedo of scattering. 

Based on a one-dimensional model of transient heat transfer, a large number of 
simulated experiments have been studied to determine which experimental conditions 
provide the best estimates of the thennal properties. 

Maximization of the optimality criterion D leads to an experiment containing one 
heating event that could be used to determine the radiation and conduction parameters. 
Strict use of the D parameter leads to experiments that were somewhat too short in the 
sense that more accurate parameter estimates could be obtained from simulated data sets 
of longer duration. By considering the conduction and radiation behavior of the model, a 



set of two shorter heating events were found that provided more accurate parameter 
estimates. 

The radiation parameter a> was easily found to a high degree of accuracy by all of 
the simulated experiments. The coupling coefficient, ‘a', was the most difficult 
parameter to determine accurately, in part because its sensitivity coefficient is similar in 
shape to that of the other conduction parameter, F. This indicates that the conduction 
parameters are close to linearly dependent. Coupling coefficient ‘a’ was an ad-hoc 
parameter added to Sullins' model in order to obtain better agreement with experimental 
data at high pressures. This suggests that there may be a different pair of conduction 
parameters that could adequately model the conduction process and provide more robust 
estimation than exhibited here. Complete results are given in Appendix C. 


Conclusions. 

Optimality criterion D is an important tool for exploring possible experimental 
conditions at little computational cost. Criterion D, however useful, is not sufficient to 
determine how well an experiment design will actually perform. Estimation of the 
parameters from simulated experimental data, or, from actual laboratory experiments, is 
required to quantify the experimental accuracy. 

Optimal experiment design involves criterion D which is computed from 
sensitivity coefficients. Sensitivity coefficients are required anyway as part of the data 
analysis to extract thermal properties from experimental data. Thus, although some time 
and effort are required for optimal experimental design, no new computational tools are 
required. 

The methods of optimal experiment design presented in this report are general and 
apply to different materials, different experimental techniques, and different experimental 
conditions. Careful use of the techniques of optimal experiment design will decrease the 
amount of time needed in the laboratory and increase the accuracy of thermal property 
estimates. 
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Abstract 

This paper is a study of optimal experiment design applied to the measure of 
thermal properties in functionally graded materials. As a first step, a material 
with linearly- varying thermal properties is analyzed, and several different tran- 
sient experimental designs are discussed. An optimality criterion, based on sen- 
sitivity coefficients, is used to identify the best experimental design. Simulated 
experimental results are analyzed to verify that the identified best experiment 
design has the smallest errors in the estimated parameters. This procedure is 
general and can be applied to design of experiments for a variety of materials. 
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Nomenclature 


bk k th parameter 

C specific heat per unit volume, [J(m 3 i^)~ 1 ] 
D optimality condition, Eq. (12) 
e slope of spatial variation, unitless 

k thermal conductivity \W{mK)~ 1 } 

L sample thickness, [m] 
n number of time steps 

p number of parameters 

qo applied heat flux, [W m“ 2 ] 

s number of sensors 

t time, [s] 

T temperature, [K] 
x spatial coordinate [m] 

X sensitivity coefficient, Eq. (9) 

X sensitivity matrix [sn x p] 

Greek 

a thermal diffusivity [m 2 s _1 ] 
e small value for finite difference 

p density [kg m“ 3 ] 

6 dimensionless time 

Subscripts 

h heater 

i index for time step 

j index for sensors 

k index for parameters 

Superscripts 

( ) spatial average quantity 

( ) + dimensionless quantity 
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Introduction 


Functionally graded materials are being studied as possible components of aero- 
space thermal protection systems. Functionally graded (FG) materials include 
composites with epoxy and metal matrices, metal foams, or any structure with 
properties designed to vary with position. In the future when FG materials 
are specified as part of a vehicle program, part of the procurement process will 
involve certification that the material meets the specifications. 

To date there has been little research on accurate thermal characterization of 
FG materials. The present research is intended to close this gap in the procure- 
ment cycle by developing accurate methods to measure the thermal properties 
of FG materials. 

A review of the pertinent literature is given next. The focus here is on 
FG materials described by macroscopic or effective properties, rather than on 
microscopic structures. Several researchers have found exact analytical so- 
lutions for thermal response by representing a FG material as composed of 
multiple layers each with different, spatially uniform, thermal properties 1,2 , 
or with exponential-function variation of thermal properties along one spatial 
direction 3,4 . Still others have used Galerkin’s method to find temperature in 
materials with arbitrary property distributions 5 . The primary motivation for 
these studies of temperature has been to determine the thermal stresses. There 
has also been some work to find the distribution of thermal properties that 
optimizes the thermal stress distribution 6,7 . 

There is only one research group that has reported experiments to measure 
thermal properties in a FG material. Makino and Noda have used transient the- 
ory applied to a FG material with exponential variation of thermal properties 8,9 . 
They have measured the single parameter that describes the thermal property 
variation in the material with a transient heating experiment. Their data anal- 
ysis combines a single temperature datum with their transient theory to provide 
a single value for the parameter. Although simple in concept, this approach is 
sensitive to measurement noise. 

The approach used in the present work is parameter estimation, a statistics- 
based method of property measurement, that has been applied to transient 
experiments for many years 10 . In this method the desired parameters are found 
by non-linear regression between the experimental data (temperatures in this 
case) and a computational model of the experiment. Parameter estimation 
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concepts have recently been applied to optimal experiment design for thermal 
characterization of uniform materials 11,12 and for materials with temperature 
varying properties 13 . 

The focus of this paper is to develop optimal experiments to thermally char- 
acterize FG materials with low thermal conductivity. To the author’s knowledge 
this is the first study of optimal experiments for thermal properties in FG ma- 
terials. 

Next a brief overview of the paper is given. In the next section a heat 
transfer model of the one-dimensional FG material is described. The sensitivity 
coefficients and sensitivity matrix are then defined, and their use in the design of 
optimal experiments is described. Several experimental designs are investigated 
for a material with linearly- varying thermal properties under specific (simu- 
lated) experimental conditions. The best experimental design and the optimal 
operating conditions are identified. The results are also verified with simulated 
experiments for estimation of thermal parameters from noise-containing data. 


Model 


In this section a one-dimensional heat conduction model of the FG material is 
discussed. This model is used to both simulate the experiments and to construct 
the sensitivity coefficient matrix. 

Consider the one-dimensional heat conduction in a slab of thickness L. The 
thermal conductivity and (volume) specific heat vary with coordinate x. The 
material properties are constant with respect to temperature (or small changes 
in temperature are assumed). The following dimensionless variables will be used 
to describe the heat conduction problem: 


i 

F< 

II 

+ 

E-i 

+ x net. 

x + ; 9 = — u 

(i) 

q 0 L/k 

L L 2 

k+ = k(x)/k- C+ = 

C{x)/C\ a — k/C 

(2) 

i r L 

— 1 f L , , 


k = — 1 k(x)dx\ 
L Jo 

C = j C(x)dx 

L Jo 

(3) 


Here k and C are the spatial average properties over the slab body, T 0 is a 
fixed temperature, and q 0 is the applied surface heat flux. The use of spatial- 
average properties to normalize the problem facilitates comparisons between 
different experiments and between different materials. 


4 


Based on the author’s previous experience with low-conductivity materials 14 , 
there are several elements that every experiment should contain: rapid heating 
on one side of the material for some period; continued data collection during a 
zero-heating period; and, a fixed temperature at the other side of the material 
(if active cooling is not practical, a large thermal mass at the non-heated face 
can be used). 

Using the above dimensionless values, this type of experiment can be simu- 
lated by solving the following equations: 


3 fV OT ’ + '|_C+ 8r+ - 

dx+ V OX+ ) 36 ’ 

0 < 

x + < 1 

(4) 

at — 0, - & + (0)— - r 

dx+ 


o < Oh 

6 > 9h 

(5) 

at x+ = 1, T + (1,0) 

= 0 


(6) 

at 0 = 0, T + (x + , 0) 

= 0 


(7) 


Heating takes place at surface x+ = 0 until time Oh after which the heating 
ends. This problem is solved by a finite difference procedure. The time deriva- 
tive is treated with the Crank- Nicholson method with uniform time steps. The 
spatial nodes are crowded towards x = 0 with a sine-squared scheme so that the 
early-time temperature can be accurately computed with a reasonable number 
of nodes. Properties k+(x) and C + (x) are evaluated in a subroutine so that 
different property distributions may be easily studied. The numerical solution 
was verified by comparison with two exact solutions: a constant-property tran- 
sient solution 15 and a steady-state solution for a material with linearly- varying 
properties 16 . These comparisons show that 40 spatial nodes are adequate and 
that the maximum timestep should be about 56 = 0.005 for 0.1% accuracy in 
the surface temperature value. 


Optimal Experiment Design 


Sensitivity coefficients and the sensitivity matrix are needed in the design 
of optimal experiments for thermal property evaluation. The sensitivity coeffi- 
cients are defined by 


Xjk(i) — bk 


dT+ 

db k 


(8) 
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which is the sensitivity for the kth parameter, the j th temperature sensor, and 
the ith time step. Parameters b may include conductivity, specific heat, density, 
etc. In this research the sensitivity matrix may encompass several heating events 
considered together as one experiment, in which case additional heating events 
are treated as additional sensors. 

The sensitivity coefficients were computed with a finite-difference procedure 
to approximate the derivative, as follows: 

v , [T${(l + e)b k )-T+(b k )] 

Xjk(i) « 6 *— ( 9 ) 


Here is the temperature at the ith timestep for the j th sensor. The value 
of e = 0.001 was found to give well-behaved values for X. 

Much can be gained from studying the sensitivity coefficients to guide the 
design of an experiment, and there are two specific requirements that the sen- 
sitivity coefficients must satisfy. First, the sensitivity coefficients should be as 
large as possible. Generally any change in the experiment that increases the size 
of the sensitivity coefficient is an improvement. Second, when two or more pa- 
rameters are to be measured in the same experiment, the sensitivity coefficients 
must be linearly independent. That is, the shape of the sensitivity coefficients 
must be different. A formal procedure to quantify these two requirements is 
given next. 

The sensitivity coefficients are assembled into a sensitivity matrix X , defined 

by 




X(l) 


X = 

X(2) 



. X W . 

X n (*) 

Xi2 (*) ••• 

Xi P (i) 

Xai(t) 

X 22 (*) ••• 

x 2p (i) 

X tl (<) 

X s2 (i) ••• 

x sp ( i) 


( 10 ) 


Optimum experiment design is based on maximization of a quantity con- 
structed from the sensitivity matrix multiplied by its transpose, given formally 


6 



( 11 ) 


by X T X. The elements of \p x p] matrix X T X are given by 

s n 

(x T x) lm = '£ / '£ x jl ( i )x jm (i) 

j— 1 i=l 


The optimality criterion selected for this study is the (normalized) determinant 
of matrix X T X , given by 





det(X T X) 


( 12 ) 


Note that the optimality criterion D is normalized by the maximum tempera- 
ture rise (squared), the number of sensors, and the number of time steps. This is 
important so that the optimality criterion D may be used to compare different 
experiments. The determinant is computed for any value of p with well-known 
matrix methods 17 . The optimality criterion insures that the sensitivity coeffi- 
cients will be large and linearly independent. 

The optimality criterion is subject to the following standard statistical as- 
sumptions: additive, uncorrelated errors with zero mean and constant variance; 
errorless independent variable; and, no prior information. Maximizing the op- 
timality criterion minimizes the hypervolume of the confidence region of the 
parameter estimates 13 . 


Experimental Designs Considered 

The focus of this research is to explore the design of experiments for ther- 
mal characterization of functionally graded materials, that is, materials with 
spatially-varying thermal properties. 

As a first step, a material with linearly varying properties was analyzed. 
Consider a one-dimensional slab body (0 < x < L) of this material. The thermal 
conductivity [W (mK)^ 1 ] and volume specific heat [J(m 3 i i C)“ 1 ] are given by 

k(x) = fc[ 1 -j- e • {x/L — 1/2)] (13) 

C{x) =C[ 1 + e • {x/L - 1/2)] (14) 

Here the parameters are k ) the spatial- average thermal conductivity, C, the 
spatial- average specific heat, and e, the dimensionless slope. The same slope is 
used for both k and C to represent the effect of density variation on thermal 
properties in a metal foam, for example. 
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Figure 1: Spatial variation of thermal properties studied. 

A finite-difference computer routine was written to compute temperature, 
sensitivity coefficients, and optimality criterion D. Three levels of spatial-property 
spatial variation were studied: e — 0.2, 1.0, and 1.8, representing property varia- 
tion of ±10%, ±50%, and ±90%, respectively, from the mean value, as shown in 
Fig. 1. Different mean values for k and C were not studied since the normalized 
results are valid for any conductivity and specific heat. 

Several combinations of simulated experimental conditions were studied, in- 
cluding the number of and location of sensors, heating on one side or the other, 
heating duration, experiment duration, and the number of parameters. 


Results for Optimality 

The results for optimality condition D are given in this section. The special 
case when the thermal properties are spatially constant is considered first for 
comparison with earlier work. There are only two parameters present, k and 
C. Consider an experiment with a single on-off heating event and with one 
temperature sensor located on the heated surface. The normalized optimality 
condition for this case is plotted versus dimensionless time in Fig. 2 for several 
different heater-off times. This figure reproduces the results of Taktak et al. 11 
to provide verification of the code used for the present research. Note in Fig. 2 
that continuous heating creates a single baseline value for each case, and then 
when heating stops the D-value jumps above this baseline by a factor of two or 
so. The optimal experiment for a uniform-property material involves a heating 
duration of Oh = 2.25 and experiment duration about 3.0 where D max = 0.02* 
In the next sections materials with linearly- varying properties will be discussed. 
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Figure 2: Optimality condition D for a uniform property material for estimat- 
ing parameters k and C for several values of heater-off time Oh . There is one 
temperature sensor at x = 0 and heating at x = 0. 


One heating event, one sensor. 

In this section a material with linear properties described by Eqs. (13) and (14) 
is studied for which the three parameters are k, C, and e. Consider a simulated 
experiment containing a single heating event and a single sensor at the heated 
surface. The other surface is maintained at ambient temperature. Consider first 
continuous heating to investigate the baseline values of D. In Fig. 3 optimality 
condition D is plotted versus time for three values of slope e for a material 
heated at x = 0 (the low-A: side). The central result shown in Fig. 3 is that 
the baseline D - values for three parameters are orders of magnitude less than for 
two parameters shown in Fig. 2. Clearly is it more difficult to estimate three 
parameters compared with two. Another observation is that the magnitude of 
D is smaller for smaller values of slope e. Thus it is more difficult to estimate 
small values of e for which the properties are nearly uniform. 

Next consider Fig. 4 in which D-values are shown for the same materials 
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at/ L 2 

Figure 3: Optimality condition D for a linearly- varying property material with 
continuous heating at re = 0 and one sensor at x = 0. Parameters are fc, C, and 


but with the experiment reversed: continuous heating at x — L (the high -k 
side); a single temperature sensor at x = L\ and, a fixed temperature at x = 0. 
For e = 0.2 the D- values are similar in size and shape to Fig. 3. For e = 1.0 
the peak D-value occurs at a slightly later time because of the greater thermal 
mass near x = L. (Recall that the time axis is normalized by the spatial- aver age 
thermal diffusivity.) Finally for e = 1.8, not only is the peak D - value delayed 
but the peak value is about 10 times greater than for heating at x = 0. At 
this point one might conclude that heating at the high-/c side is best, at least 
for large-e materials. However larger £>-values than these are presented in the 
following sections by the use of interior temperature sensors and by combining 
two heating events in a single experiment. 

One heating event, two sensors 

In this section simulated experiments were analyzed with two sensors, one at the 
heated surface and one sensor inside the sample. Generally a second temperature 
sensor located in the range 0.2 < x/L < 0.4 gave the largest D- values. If the 
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at/L 2 

Figure 4: Optimality condition D for a linearly- varying property material with 
continuous heating at x = L and one sensor at x = L. Parameters are k , C, 
and e. 

second sensor is too close to the surface sensor it offers little new information, 
and if the second sensor is too close to x = L its response will be limited by the 
fixed- temperature boundary there. 

In all of the cases reported here the second sensor is located at x^ — 0.25L. 
Results for e = 0.2 are typical and are presented in Fig. 5. Note that the 
baseline values for D for two sensors, with the heater continuously energized, 
are about 60 times larger than for one sensor as shown in Fig. 3 (for case 
e = 0.2). That is, use of an additional sensor inside the body greatly improves 
the experiment. In Fig. 5 four cases are also presented for which the heater is 
shut off before 0 = 4.0. The best experiment is that for which the heater is shut 
off at 0^ — 2.3 and the experiment continues until 0 — 2.8. The same trend of 
improvement is present for other values of e. Results for e = 1.0 and e = 1.8 
for two sensors and one heating event are listed as experiment 1 in Table 1. 
The table shows that the optimal heating times and experiment durations are 
shorter for larger values of e. For all e- values the choice of heater shut-off time 
causes only small changes in the maximum D- value, that is, the peak D~ value 
is insensitive to heating duration when an interior sensor is used. 
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Table 1. Maximum value of optimality condition D, and the conditions under 
which it occurs, for three experiments: (1) one heating event and two sensors 
at x/L = 0, 0.25; (2) two heating events, from each side, with one sensor on the 
heated side; and (3) spatially uniform properties (for comparison), one heating 
event, one sensor at x = 0. 


Experiment 

slope 

time 

time at 


(parameters) 

e 

Oh 

Dmax 

value 

1. (k, C, e) 

0.2 

2.3 

2.8 

1.7(10~ 6 ) 


1.0 

1.5 

1.8 

8.7(10 -5 ) 


1.8 

1.5 

1.7 

1.2(10 -3 ) 

2. (Jc, C,e) 

0.2 

1.3 

1.9 

9.5(10 -6 ) 


1.0 

1.7 

2.2 

24. (10 -5 ) 


1.8 

1.3 

1.5 

2.2(10 -3 ) 

3. (k, C) 

- 

2.25 

3.0 

0.02 


Two heating events, one sensor each. 

In this section results are presented for an experiment composed of data from 
two heating events, each involving one surface-mounted temperature sensor. In 
one heating event the sample is heated at x = 0 and the sensor is located at 
x — 0. In the other heating event the sample is heated at x = L and the 
sensor is located at x — L. In the laboratory the second heating event could be 
accomplished with the same heater and sensor by reversing the sample. This 
experimental design takes advantage of the different conductivity values on each 
side of the body and the fact that surface-mounted sensors are simpler to install 
than interior sensors. 

For simplicity both heating events involve the same heating duration 0 h 
and the same data duration. In Fig. 6 D-values are plotted versus time for 
four different heater-off times, again for e = 0.2. Each curve represents the 
combination of data from two heating events into a single experiment. Once 
again the familiar shape occurs with a baseline value for continuous heating, 
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Figure 5: Optimality condition D for an experiment with two temperature 
sensors at x/L ~ 0 and 0.25 and heating at x = 0 for several cases with different 
heating duration. Parameters are fc, (7, and e. 

with heater-off cases providing an additional boost to the maximum D value. 
The distinguishing feature, once again, is the magnitude of D max compared to 
earlier experiment designs. This case with two heating events gives a D max 
nearly 6 times higher than for the interior sensor case (Fig. 5). The best 
experiment from the e = 0.2 results shown in Fig. 5 is for heating off at Oh = 1.3 
which provides D m ax = 9.5(10 -6 ) at experiment duration 0 = 1.9. 

A summary of results for e = 1.0 and e = 1.8 for two heating events are 
listed as experiment 2 in Table 1. The same general trends are exhibited for 
these e-values, however the amount of improvement in Dmax is less for higher 
e- values when comparing two heating events with an interior sensor. 


Simulated Experiments 

The purpose of optimal experimental design is to provide meaningful assis- 
tance in thermal characterization and in analyzing experimental data. In this 
section simulated thermal-characterization experiments are carried out for two 
experimental designs. 


13 



Figure 6: Optimality condition D for two heating events , one at x = 0 and one 
at x = L, combined into one experiment. The temperature sensor is located at 
the heated surface. Parameters are k, C, and e. 

Simulated experiments are carried out by adding errors to exact temperature 
values (computed from the model), and analyzing this error-containing data to 
estimate parameters. The added errors are normally distributed with zero mean 
and an adjustable variance. The error variance was set to either 1% or 5% of the 
maximum temperature values. The added error values are found with a com- 
puter routine that requires a seed number, and different seed numbers can be 
used to produce different sequences of error values (these are sometimes called 
pseudo-random numbers). The simulated data is analyzed with a Marquardt 
regression scheme which systematically compares the simulated data with val- 
ues computed from the model based on iteratively improved guesses for the 
parameters. Iteration ceases when the improvements in the parameters become 
small. 

Figure 7 shows a set of simulated data with added error variance 5%. The 
regression fit for this data is also shown in the figure. The 5% variance in the 
added error is much larger than would be tolerated in reasonable experimental 
practice, but it is useful as a test of the estimation scheme. In this case the 
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Figure 7: Simulated data (with noise) and regression fit for experiment 1 (one 
heat event, two sensors) for e = 1.8 and error variance 5%. 


estimation scheme converges without incident. 

The regression scheme was applied to a variety of simulated experiments, 
and it was found that the parameter estimates varied somewhat for different 
error sequences (different seed numbers in the random- number generator), even 
though the variance of the errors was identical. This suggested that a single 
simulated experiment may be misleading as to the precision of the estimates. 
One estimate might by chance be particularly close to the actual value, and 
another estimate might be particularly far from its actual value. To deal with 
this uncertainty, each simulated experiment was repeated ten times, and an 
ensemble average error was computed. For each repetition n, the error for 
parameter b t is given by: 


erri(n) = 


b% b% 

hi 


x 100% 


(15) 


Here the exact value is 6* and the estimated value is b l . The ensemble average 
error for parameter bi is given by 

1 10 

erf< = 53 |eTTi(n)| (16) 


n—1 
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Note that the absolute value of each error is used in the ensemble average. The 
purpose here is to compare the errors produced by different experiments, and 
the absolute value reduces the variability in the results without biasing the error 
values towards zero. In contrast, in analyzing lab data from several identical 
experiments, a simple average of the parameter estimates (no absolute value) 
would be appropriate to find the most precise estimates. 

Table 2 shows the summary of results of the ensemble average error found 
from ten repetitions of each simulated experiment. Two experiment designs are 
compared in the table, and the operating conditions for each experiment axe 
taken from the optimal conditions given earlier. The values in Table 2 show 
that experiment 2 tends to provide smaller error in the estimates for the high- 
noise data, and when the material distribution slope, e, is small. Under these 
difficult conditions, the conclusion is that experiment 2 is best for estimating 
parameters. For less difficult conditions, for example with low-noise data and 
for e large, experiments 1 and 2 provide comparable-size errors in the estimated 
parameters. 

In all cases listed in Table 2 the simulated experiments were carried out with 
30 simulated data points (extracted from the many time steps required in the 
model calculation), and the initial guesses for the parameters in the regression 
scheme were taken to be 0.9 times the correct parameter values. To explore the 
effect of these arbitrary choices on the parameter estimates, additional simulated 
experiments were carried out. Additional cases included: 60 simulated data 
points; initial guesses of 0.5 times the correct parameter values; and, initial 
guesses of 1.5 times the correct parameter values. In all these additional cases 
the results were comparable to those in Table 2. 

Summary and Conclusions 

In this paper optimal experiments are sought for the measurement of thermal 
properties in FG materials with spatially-varying thermal properties. As a first 
step, one- dimensional transient experiments were studied for a material with 
linearly-varying thermal properties described by parameters C, fc, and slope e. 
Variation of properties with temperature was not treated. 
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Table 2. Percent error in parameter estimates from simulated experiments. Re- 
ported values are ensemble averages over 10 repetitions of the data analysis for 
each experiment design. 



percent 


ensemble average error 

Experiment 

variance 

slope 

in parameters, percent 

Design 

in data 

e 

k 

C 

e 

1. One heat event, 

1 

0.2 

0.2435 

2.1681 

10.9687 

two sensors at 

1 

1.0 

0.3100 

1.9259 

1.2855 

x/L = 0, 0.25. 

1 

1.8 

0.2936 

1.5820 

0.1461 


5 

0.2 

0.9914 

10.6676 

59.2553 


5 

1.0 

2.2731 

13.4547 

18.2089 


5 

1.8 

1.9463 

4.7742 

0.9925 

2. Two heat events, 

1 

0.2 

0.3249 

1.0254 

7.4094 

from each side, 

1 

1.0 

0.3706 

0.8316 

1.3456 

one sensor on 

1 

1.8 

0.8085 

0.5851 

0.3114 

heated side. 







5 

0.2 

1.5017 

6.8977 

23.6133 


5 

1.0 

1.9755 

4.0834 

6.3749 


5 

1.8 

3.1122 

1.7011 

1.4480 


Several transient experimental designs were considered, including single or 
multiple heating events and including surface and/or interior temperature sen- 
sors. An optimality condition based on sensitivity coefficients was used to find 
the best operating conditions for each experiment, and to find the best exper- 
imental design among those studied. The best experiment has the smallest 
hypervolume of the confidence region of the parameter estimates. 

The results of the optimality study show that it is more difficult to obtain 
accurate values of slope e when e is small compared to when e is large; that is, 
it is easier to “see” large spatial variations in properties. The best experimental 
design involves analysis of combined data from two separate heating events, 
one with heating on one side of the body and one with heating on the other 
side, each time with the sensor located on the heated side and the unheated 
side maintained at a fixed temperature. This conclusion is also supported by 
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a series of simulated experiments, carried out with regression analysis of error- 
containing temperature values. 

The optimal operating conditions for the experiment are somewhat depen- 
dent on the property slope e. For the specific case e = 0.2, the optimal conditions 
are heating duration of Oh = 1.3 and with data recorded until 0 — 1.9, where 6 
is dimensionless time. 

The method of optimal experiment design discussed in this paper has been 
demonstrated for a class of functionally- graded materials, however the approach 
is completely general and applies to any material. Work in progress includes op- 
timal experiment design for thermal characterization of porous materials under 
conditions where both conduction and radiation heat transfer are present. 
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ANALYSIS OF PHOTOTHERMAL CHARACTERIZATION 

OF LAYERED MATERIAL S- 

DESIGN OF OPTIMAL EXPERIMENTS 1 

Kevin D. Cole 2 


Abstract 

In this paper numerical calculations are presented for the steady-periodic tempera- 
ture in layered materials and functionally-graded materials to simulate photothermal 
methods for the measurement of thermal properties. No laboratory experiments were 
performed. The temperature is found from a new Green’s function formulation which 
is particularly well-suited to machine calculation. The simulation method is verified 
by comparison with literature data for a layered material. The method is applied to a 
class of two-component functionally-graded materials and results for temperature and 
sensitivity coefficients are presented. An optimality criterion, based on the sensitivity 
coefficients, is used for choosing what experimental conditions will be needed for pho- 
tothermal measurements to determine the spatial distribution of thermal properties. 
This method for optimal experiment design is completely general and may be applied 
to any photothermal technique and to any functionally-graded material. 

KEY WORDS: functionally graded material; Green’s functions; optimal experiment; 
photothermal; thermal properties 


1 Introduction 

Functionally-graded (FG) materials are being studied as possible components of aero- 
space thermal protection systems. These materials include composites with epoxy and 
metal matrices, metal foams, or any structure with properties designed to vary with 
position. In the future when FG materials are specified as part of a vehicle program, 
part of the procurement process will involve certification that the material meets the 
specifications. 

To date there has been little research on accurate thermal characterization of FG 
materials. The present research is intended to close this gap in the procurement cycle 
by investigating photothermal methods for non-destructive and accurate measurement 
of thermal properties in FG materials. In this paper only numerical simulations are 
presented and no laboratory experiments were performed. 

1 Paper presented at the Fifteenth Symposium on Thermophysical Properties, June 22-27, 2003, 
Boulder, Colorado, U.S.A. 

2 Mechanical Engineering Dept., University of Nebraska-Lincoln, Lincoln, NE 68588-0656 
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A review of the pertinent literature is given next in three areas: computer simu- 
lation of FG materials; heat transfer theory for photothermal methods; and, optimal 
experiment design. 

Several researchers have found analytical solutions for the thermal response by 
representing a FG material as composed of multiple layers each with different, spa- 
tially uniform, thermal properties [1,2,3]. Other studies investigated an exponential- 
function variation of thermal properties along one spatial direction [4,5]. Another used 
Galerkin’s method to find temperature in materials with arbitrary property distribu- 
tions [6]. The primary motivation for these studies has been to determine temperature 
for the purpose of finding thermal stresses, or to find the distribution of thermal prop- 
erties that optimizes the thermal stresses [7,8]. 

One research group has reported transient-heating experiments to measure thermal 
properties in a FG material [9,10]. This group studied a FG material containing an 
exponentially-varying spatial distribution of thermal properties. Their data analysis 
combines a single temperature datum with their transient theory to provide a sin- 
gle value for the parameter describing the spatial distribution of thermal properties. 
Although simple in concept, this approach is sensitive to measurement noise. 

In the area of photothermal measurements, there are several pertinent publica- 
tions. A diverse collection of thermal-wave Green’s functions and temperature solu- 
tions has been published recently in book form [11]. Primarily homogeneous materials 
are treated, and layered materials are included by defining a global Green’s function 
that embodies the effects of several layers in the material. Because the complexity of 
the layered-body Green’s function increases rapidly as layers are added, no more than 
3 layers are discussed. 

Theory for many-layered bodies to laser heating has been studied previously by the 
author [12,13]. The volumetric heating is treated exactly from the optical absorption 
properties of all layers. The multi-layer body is treated efficiently by the use of local 
Green’s functions which are found first in the time domain and are then transformed 
into the frequency domain. Each layer is linked to adjacent layers with appropriate 
interface conditions. 

Theory for the photoacoustic response of a layered solid has been recently studied 
for which the temperature in each layer is linked with adjacent layers by interface 
conditions [14]. The optical absorption in each solid layer is described by an exponen- 
tial distribution and an absorption coefficient. The photoacoustic response is found 
from both thermal effects and mechanical effects in the gas, however thermal effects 
predominate for solid materials. The method is used for analysis of experimental data 
in materials with 2 and 3 layers. 

In the area of optimal experiment design, parameter estimation has been used for 
obtaining thermal properties from transient experiments for many years [15]. In these 
methods the desired parameters are found by non-linear regression between the experi- 
mental data (temperatures in this case) and a computational model of the experiment. 
Parameter est im ation concepts have recently been applied to optimal experiment de- 
sign for thermal characterization of uniform materials [16] and for materials with 


2 



temperature- varying properties [17]. The author has previously studied optimal ex- 
periment design for low-conductivity FG materials [18]. The simulated experiments 
involved time-series data collected from one or more temperature sensors, and the 
data analysis is carried out in the time domain. An optimality criterion was used to 
find the best experimental conditions for simultaneous estimation of several thermal 
properties. Results of simulations show that for FG materials with spatially-varying 
conductivity, it is better to heat the sample from the low-conductivity side. 

In the present paper, FG materials are simulated with a large number of intercon- 
nected layers. The heat transfer theory draws upon the author’s previous work with 
Green’s functions, but here the Green’s functions are given directly in the frequency 
domain in the form of algebraic expressions, not infinite-series expressions, that are 
numerically well behaved under all conditions. Likewise the temperature expressions 
found from these Green’s functions are numerically well behaved. The Green’s func- 
tions are given for a variety of boundary conditions; previously only specified-flux 
boundaries were treated. To the author’s knowledge this paper describes the first 
application of optimal experiment design methods to frequency-domain analysis ap- 
propriate for photothermal experiments. This paper is divided into several sections, 
as follows: the temperature in one layer; the Green’s functions; the temperature in a 
multi-layer material; the design of optimal experiments; results for a layered material; 
results for a FG material; and, a brief summary. 


2 Temperature in one layer 

Since the photothermal applications of interest involve periodic heating by a laser, the 
solution is sought in Fourier-transform space, and the solution is interpreted as the 
steady-periodic response at a single frequency ui. For a discussion of this point see 
[19]. Consider the heat conduction equation in Fourier transform space in one layer: 

§--f r = 0<X<L (1) 

dT 

ki~ VhiT — fi{u)-, at boundaries i = 1,2 (2) 

orii 

Here T is Fourier-space temperature (K s), a is thermal diffusivity (m 2 /s), k is thermal 
conductivity (W/m/K), g is volume heating (W s/m 3 ) deposited by a laser, and /* is a 
specified boundary condition. Index i = 1,2 represents the boundaries at the limiting 
values of coordinate x. The boundary condition may be one of three types at each 
boundary: for type 1 fo is a specified temperature (k{ = 0 and hi = 1); for type 2 /, is 
a specified heat flux (hi — 0); and, type 3 represents a convection condition where h t 
is a constant- with-time heat transfer coefficient (or contact conductance). 

The temperature will be found with the Fourier-space Green’s function, defined by 
the following equations: 

~ ° 2 G = -~&(x-x') (3) 

ox z a 
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(4) 


ki~ 1- hiG = 0; i — 1,2 

OTli 

Here a 2 = ju/a and 6(x-x') is the Dirac delta function. The coefficient 1/a preceding 
the delta function in Eq. (3) provides the frequency-domain Green’s function with 
units of seconds/meters. This is consistent with our earlier work with time-domain 
Green’s functions. 

Assume for the moment that the Green’s function G is known, then the steady- 
periodic temperature is given by the following integral equation (see [20], p. 40-43): 


T{x,uj) 


J g(x', cj)G(x, x',u)dx' (for volume heating) 


+afi(u>) x 


dG/dn'(x,Xi,oj) 

\G{x,Xi,u)) 


(type 1 only) 
(type 2 or 3) 


i = 1,2 


( 5 ) 


3 Green’s Function 


The Green’s function (GF) that satisfies Eqs. (3) and (4) is given by 


G(x,x',uj) = 


+ 


S 2 (ST e~ <r ( 2L - |x - x 'l ) + 5 1 + e~ <T(2L_x_a: ' ) ) 
2aa{S/St - SfS^eF^ 1 ) 
Sj(5fe-^ | *"* ,|) + Si e -*(*+*'>) 

2aa{S± S} - 5f5 2 -e-^) 


( 6 ) 


where the subscripts 1 and 2 represent the two boundaries at the smallest and largest 
x-values, respectively. Coefficients S m and depend on the boundary conditions on 
side M and are given by 


‘S’m 




1 

kc T + hf4 


if side M is type 0, type 1, or type 2 
if side M is type 3 


0 if side M is type 0 

—1 if side M is type 1 

1 if side M is type 2 

k ka — h M if side M is type 3 


A boundary of type 0 designates a far-away boundary, as in a semi-infinite body. The 
derivation of the Fourier-space GF in Eq. (6) parallels that for steady-state GF given 
elsewhere [21]; however in the present work a is complex. 

This form of the GF is particularly well-behaved for machine computation, and 
most importantly, the temperature expressions based on these GF are similarly well- 
behaved for any layer thickness and for any frequency. This is a key contribution of this 
paper, in sharp contrast with previously reported difficulties in evaluating numerical 
values from exact solutions. For example, numerical overflow can occur with other 
formulations in thermally thick layers [14]. In a time-domain study, only short-time 
results were included due to numerical difficulties associated with longer times [3]. 
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The GF expression given in Eq. (6) covers a number of boundary condition com- 
binations, and a numbering system is used to distinguish among them. Designation 
XIJ is used to identify the GF for heat transfer in a layer with boundary condition 
of type I = 1, 2, or 3 at x = 0 and with boundary condition of type J = 1, 2, or 3 
at x = L. For example, designation X12 represents the GF with type 1 boundary at 
x = 0 and type 2 boundary at x — L. Designation XIO is used to identify the GF for 
a semi-infinite region with a boundary of type I = 1, 2, or 3 at a; = 0. 

4 Temperature in Layered Materials 

In this section the temperature caused by absorption of laser energy will be found 
in a domain consisting of non-absorbing air, N solid layers, and a substrate. At 
the interfaces, let q nm represent the heat flux leaving layer n and entering layer m. 
Applying Eq. (5), the interface temperature in the air is: 

T o (0, w) = — ^Go(0, 0, u)q w (7) 

Kq 

In layer v, i — 1,2, N : the interface temperatures axe: 

Ti(0,u>) = -j^Gi(0,0,u)qi-i t i + -j^Gi(0,Li,u})q i+ i )i + Bi(0) (8) 

r V| rCj 

Ti(Li,u ;) = -£-Gi(Li,0,u)qi-ij + ^-Gi(Li,Li,u)qi + i t i + Bi(Li) ( 9 ) 

/Cj Ki 

In the substrate the temperature at the interface is: 

Tv +1 (0,u;)= Gn+i (0, 0, ^)gw,N + i + -Efyr+i (0) (10) 

«w+i 

In the above expression, symbol Bi has been used for the volume-heating integral term 
from Eq. (5), specifically, 

Bi{x) = -r [ g{x',u) Gi(x,x',u}) dx' (11) 

Ki Jx* 

Here g is the laser energy absorbed in the layer per unit volume; this can be determined 
without approximation from the optical properties of the layers [12], 

In the above temperature expressions, all of the interface heat fluxes are initially 
unknown. The heat flux leaving one layer enters the adjacent layer, = —q^i-x 

and the temperature difference between adjacent layers is related to heat flux through 
a contact resistance at each interface: 

Qi-i,i Ri = Ti(0, cj) — Ti-i(Li-i, u); i — 1, 2, . . . , N + 1 (12) 

The contact resistance Ri describes the size of the temperature jump across the inter- 
face. Next Eqs. (7-10) are combined with Eq. (12) to eliminate temperature. The 


result is a set of N + 1 linear algebraic equations for the unknown heat fluxes, which 
may be stated in matrix form: 

Wo U i R\ — Vi 0 ... 0 

—Vi Ui + U 2 + R 2 — V 2 ••• 0 

0 —V 2 U 2 + U 3 4- ••• 0 

-V N 

0 0 ... —Vff Un + W/v+i + Rn+1 



010 


' ^(0) 


921 


B 2 ( 0 ) - Bi(L\) 

X 

932 

— 

B s (0) — B2(L 2 ) 


. Qn+i,n . 


. Bn + i(0) - Bh(Ln) 


Symbols Wi, Uj, and V* used in the above expression are given below: 

w, = 2t G( (0,0,w) 

U, = %G l (0,a,u)=%G i (L l ,L„u) 

V = %-G i (0,L i ,u>) = ^G i (L i ,0,u) (14) 

For any multilayered system, it is now possible to calculate the heat fluxes (q t j) through 
all interfaces in the system. The above result is exact, and Cramer’s rule may be used 
to solve for the q's for a sample composed of one or two layers. For a sample with two 
or more layers, a numerical matrix solution is best. Once the heat fluxes are found, 
the temperature at any interface is given by Eq. (8 - 10), or the temperature within 
any layer may be found with Eq. (5). 

Several different GF may be used in the above matrix equation. The non-absorbing 
gas (region 0) is a semi-infinite region so the GF needed is number X20. For layers 
i = 1, 2, . . ., N the GF needed are type X22 (specified heat flux). The GF for 
the substrate depends on the heat transfer environment there. For example, a thick 
substrate could be described by GF number X20, or, a substrate in imperfect contact 
with a cold plate could be described by GF number X23. 

5 Optimal Experiment Design 

Sensitivity coefficients are central to the design of optimal experiments for thermal 
property evaluation. The sensitivity coefficients are defined by 

Xrt 0 = b (15) 

which is the sensitivity for the kih parameter, the jth temperature sensor, and the «th 
frequency. Parameters 6* may include conductivity, specific heat, density, etc. In this 
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research the sensitivity matrix has been computed from the real-valued amplitude and 
real- valued phase of the complex temperature at r frequencies, for which the distinct 
amplitude and phase values are treated as 2 r measurements from each sensor. 

The sensitivity coefficients were computed with a finite-difference procedure to 
approximate the derivative, as follows: 


eb k 


(16) 


Here Ty is the temperature at the ith frequency for the jth sensor. The value of 
e = 0.001 was found to give well-behaved values for X. 

There are two specific requirements that the sensitivity coefficients must satisfy. 
First, the sensitivity coefficients should be as large as possible. Second, when two or 
more parameters are to be measured in the same experiment, the sensitivity coefficients 
must be linearly independent. A formal procedure to quantify these two requirements 
may be constructed if the sensitivity coefficients are assembled into a sensitivity matrix 
X, which is then multiplied by its transpose, given formally by X T X, of size [p x p]. 
The optimality criterion is the (normalized) determinant of matrix X T X, given by 
[16] 

D = 77(t^ det(xTx> (n) 

Note that optimality criterion D is normalized by the maximum temperature rise 
(squared), the number of sensors s, and the number of frequencies r. This is important 
so that D may be used to compare different experiments. When D is large, the 
sensitivity coefficients will be large and linearly independent [17]. 


6 Results for a Layered Material 

In this section the techniques of simulation and optimal experiment design are applied 
to a SiC >2 film on a Si substrate. This material has been studied elsewhere by a pho- 
toacoustic technique over a frequency range from 2 to 20 kHz [14]. An opaque coating 
of Ni of 20 nm thickness is added to the sample to improve the optical absorption. 
Published results for the thermal conductivity of the SiC >2 film and the contact resis- 
tance between the film and the Si substrate are: k = 1.52 W/m/K and R < 10 -9 K 
m 2 /W. 

Using the present methods, the computed phase of the temperature agreed with 
published values, thus verifying the present approach. The sensitivity coefficients for 
the phase of the surface temperature to variations in the conductivity of the the SiC >2 
layer, k , and the contact resistance with the Si substrate, R, are plotted in Fig. 1. 
The property values used are given in Table 1. The sensitivity to conductivity k has 
a maximum (negative) value at about 6000 Hz. The shape of the R-sensitivity curves 
and the location of the largest value depends strongly on the R value. For R = 10 -8 
and R = 10“ 7 , the sensitivity curves are relatively small, with a broad peak near 3 
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Figure 1: Phase sensitivity to thermal conductivity and to layer-substrate contact 
resistance for a layer of SiC >2 on silicon at various frequencies. 

kHz. For R = 10 -6 however, the sensitivity values are large and positive at small 
frequencies and the curve slopes down to negative values as frequency increases. Some 
conclusions for R < 10 -7 are that the sensitivity coefficient is small, but the frequency 
range considered (0.5 - 20 kHz) captures the peak sensitivity. For R = 10 -6 the largest 
sensitivity may lie outside this frequency range. 

Table 1. Properties used in calculations of a three-layer solid which is 
heated at the Ni surface and exposed to a layer of air on either side. 


layer 

i 

d (fi m) 

k (W/m/K) 

a (m 2 /s) 

R (K m 2 /W) 

air 

0 

- 

0.0263 

2.25E-05 

- 

Ni 

1 

0.07 

80 

1.98E-05 

0 . 

Si0 2 

2 

0.485 

1.52 

9.09E-07 

0 . 

Si 

3 

382 

151 

9.09E-04 

varies 

air 

4 

- 

0.0263 

2.25E-05 

- 


Next the sensitivities for k and R will be examined together. Optimality criterion 
D provides a numerical measure of the extent to which the sensitivity coefficients are 
both large and linearly independent, and at which frequencies this occurs. Figure 
2 shows the optimality criterion for the layered material for both conductivity and 
contact resistance considered together. The highest curve is for k and R = 10 -6 , and 
the peak for this curve is about 4 kHz. For the R = 10 -7 and R = 10 -8 curves the 
peak occurs around 8 kHz. Figure 2 indicates that the best experiment to measure 
both k and R when R = 10~ 6 includes data at 4 kHz. For smaller i?-values the most 
important frequency is 8 kHz, however because D is smaller the analysis of the data 
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Figure 2: Optimality criterion D for simultaneous estimation of thermal conductivity 
and contact resistance, versus frequency, for a layer of Si 02 on silicon. 

may be more difficult. It is important to note that each point on Fig. 2 represents 
a value for D computed from a range of data extending from 20 kHz down to that 
point; data is added from high frequency to low frequency. To repeat, the D-values 
indicate what frequency range is needed for optimal estimation of both k and R from 
experimental data. 


7 Results for a Functionally Graded Material 

In this section the methods for experimental design are applied to a two-phase ce- 
ramic/ceramic material with graded volume fraction of the components. The volume 
fraction profile is assumed to have the form: 

VU*) = 1 - ( x/LY ; V 2 (x) = ( x/Lf (18) 

where Vi and V 2 are the volume fraction of the components. At location x = 0 the 
material is pure component 1 and at x = L the material is pure component 2. The 
particular material considered is composed of TiC and SiC and the thermal properties 
are given in Table 2. 

Table 2. Material properties of TiC and SiC. 



material 

k (W/m/K) 

9 (kg/m 3 ) 

c (J/kg/K) 

1 

TiC 

20 

4900 

700 

2 

SiC 

60 

3200 

1000 
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Figure 3: Thermal conductivity distribution in a functionally-graded material for sev- 
eral values of distribution parameter p. 


The thermal conductivity of the material is given by 


k(x) = ki 


. V2(x)(k 2 -k l ) ' 

ki + ih-k^V l (x)/3 


(19) 


where subscripts 1 and 2 stand for the properties of components 1 and 2, respectively. 
Figure 3 shows the spatial distribution of conductivity for the TiC/SiC material for 
several values of exponent p. The mass density and specific heat are determined by 
the rule of mixtures: 


p{x) = V\{x)p\ -I- V 2 (x)p 2 (20) 

c{x) = Vi(*)c! + V 2 {x)c 2 (21) 

The thermal-stress behavior of this material has previously been studied [3]. 

In photothermal methods, the sample is heated by a periodically modulated laser 
beam, and the surface temperature (or a subsequent acoustic signal) is measured 
at the modulation frequency. The temperatures reported below are computed with 
the layered GF method with 50 layers used to simulate the spatial variation in the 
sample. Both surfaces of the sample are exposed to air. The surface temperature is 
shown in Fig. 4a (amplitude) and Fig. 4b (phase) versus dimensionless frequency. 
The frequency is normalized as /* = fL 2 /a av where L is the material thickness and 
ot a v = (<*1 + oc 2 )/ 2 , the average of the component values. In Fig. 4 the sample is 
heated at x = 0 and the temperature is reported at x = 0, the low-fc side. This 
heating condition provides slightly higher temperature response than for heating on 
the high-A; side (at x — L). Note that in Fig. 4a the amplitude curves are monotonic 
with few distinguishing features. In contrast the phase has a distinctive maximum 
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Figure 4: Surface temperature (a) amplitude, and (b) phase, versus dimensionless 
frequency for a functionally-graded material heated periodically at the surface. 


for each p value, which supports the experimental observation that phase is more 
important than amplitude for photothermal measurement of thermal properties. 

Figure 5a shows the sensitivity of the phase of temperature to k 2 , the thermal 
conductivity of component 2, for three values of spatial distribution parameter p. Note 
that the sensitivity is largest for p = 1, the linear ^-distribution, and is small for other 
distributions. In Fig. 5a the peak sensitivities in a range of dimensionless frequencies 
below /* = 1. This range of frequencies represents thermal waves (generated by 
the periodic heating) that penetrate all the way through the sample thickness. The 
sensitivity to k\, the component 1 thermal conductivity, is not shown because it is 
similar in size and shape to the k 2 sensitivities (but values are negative). 

Figure 5b shows the sensitivity of the phase of the temperature to spatial distri- 
bution parameter p. The curve for p = 1 has a large positive peak at about /* = 0.2 
and a negative peak at about /* = 1.5. The largest sensitivities are again located in 
the range /* < 1. 

The sensitivity to a single parameter is useful for determining a single parameter 
from an experiment. When two or more parameters are to be determined, optimality 
criterion D is instructive. Figure 6a shows optimality criterion D for both conductivity 
k 2 and exponent p. Both amplitude and phase information was used to compute these 
values. The largest curve is for p = 1 which indicates that this spatial distribution 
will provide better estimates of thermal properties k 2 and p from an experiment than 
for other p-values. The location of the peak for each curve indicates which frequencies 
should be included in an experiment. Since the D-values shown were computed from 
high-to-low frequencies, the peak indicates the lowest frequency of data that is needed 
for optimal data analysis. 






Figure 5: Sensitivity of the phase of temperature to (a) conductivity k 2 , and (b) 
exponent p, for the functionally-graded material heated periodically at the surface. 


Optimality criterion D was also investigated for other combinations of parameters. 
The D-v alues computed for p and k\ together, not shown, were identical to Fig. 6a, 
which is expected from the similarity in the shape of the sensitivities to k\ and k 2 
mentioned earlier. The O- values computed for k\ and k 2 together are shown in Fig. 
6b. The important feature of this figure is the vanishingly small values, five orders of 
magnitude smaller than Fig. 6a values, which indicates that k\ and k 2 should not be 
sought simultaneously. As mentioned before the sensitivity coefficients for ki and k 2 
have a similar shape, so the fact that their combined D - value is near zero reinforces the 
idea that the sensitivities for k\ and k 2 are not linearly independent. Finally, D-values 
were also computed for p, k\, and k 2 considered simultaneously (not shown), and these 
values were predictably near zero because of the dependence of ki and k 2 . 

8 Summary 

This paper investigates photothermal methods for thermal characterization of func- 
tionally graded materials through numerical modeling and experiment design. No 
laboratory experiments are reported. A new formulation is given for the temperature 
response of layered materials to periodic heating, based on the method of Green’s 
functions, which is numerically better behaved than previous work. The method has 
been applied to a SiC >2 layer on silicon and compared to literature values to validate 
the method. Optimality criterion D indicates which frequency range of experimental 
data should provide the best possible estimates of the layer conductivity and contact 
conductance. 

The new methods have also been applied to a two-component functionally-graded 
material with a power-law distribution of thermal properties. The largest temperature 
response is found by heating the sample on the low-k side, and the phase of the tem- 
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Figure 6: Optimality criterion D for simultaneous estimation of (a) conductivity k 2 
distribution parameter p, and (b) conductivities k\ and k 2 , for the functionally-graded 
material. 


perature is particularly important for estimation of thermal properties. Component 
conductivities k\ and k 2 have similar-shaped sensitivity coefficients, and consequently 
both cannot be estimated simultaneously from experimental data. The most impor- 
tant parameter is p which describes the spatial distribution of thermal properties in 
the functionally-graded material. Values for optimality criterion D indicate that val- 
ues for p may be found simultaneously with one of the conductivities, but not both. 
Dimensionless frequencies in the range /* < 1 are important for measurement of 
spatial-distribution parameter p. The magnitude of the optimality criterion D also 
suggests that it will be easier to estimate parameters for p rj 1 (near-linear spatial 
variation) compared to other values of p. 
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Abstract 

Metallic foams are being investigated for possible use in the thermal protection systems of reusable 
launch vehicles. As a result, the performance of these materials needs to be characterized over a wide 
range of temperatures and pressures. In this paper a radiation/conduction model is presented for heat 
transfer in metallic foams. Candidates for the optimal transient experiment to determine the intrinsic 
properties of the model are found by two methods. First, an optimality criterion is used to find an 
experiment to find all of the parameters using one heating event. Second, a pair of heating events is 
used to determine the parameters in which one heating event is optimal for finding the parameters 
related to conduction, while the other heating event is optimal for finding the parameters associated 
with radiation. Simulated data containing random noise was analyzed to determine the parameters 
using both methods. In all cases the parameter estimates could be improved by analyzing a larger 
data record than suggested by the optimality criterion. 
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Nomenclature 


a — Coupling coefficient, m K /W 3 
b = Heater thickness, m 
c — Weighted specific heat, J/(kg K) 
c g = Specific heat of gas, J/(kg K) 
c/t = Specific heat of heater, J/(kg K) 
c m = Specific heat of metal, J/(kg K) 

Dj = Fiber diameter, m 
d g = Gas collision diameter, m 

e = eO 4- ep * T, Specific extinction coefficient, m 2 /kg 
F = Efficiency factor for conduction 
K b = Boltzmann constant 
Kn — Knudsen number 
k — Thermal conductivity, W/(m K) 
k g = Gas thermal conductivity, W/(m K) 
k* — Gas thermal conductivity at 1 atm, W/(m K) 
km = Thermal conductivity of metal, W/(m K) 

6 — Characteristic length, m 

P — Pressure, Pa 

Pr = Prandtl number 
T = Temperature, K 
T s = Sampled temperature, K 
t — Time, s 
x — Spatial coordinate, m 
a — Thermal accommodation coefficient 
(3 — Specific extinction coefficient, W/(m 3 K 4 ) 

e — Porosity 

6\ — Emittance of the septum plate 

£2 — Emittance of cold plate 

7 — Specific heat ratio 

A m = Molecular mean free path, m 
q c = Conduction heat flux, W /m 2 
q r — Radiation heat flux, W/m 2 
qr = Total heat flux, W /m 2 
p g = Density of the gas, kg/m 3 
Ph ~ Density of the heater, kg/m 3 
p m = Density of metal, kg/m 3 
Pf = Density of the foam, kg/m 3 
a — Stefan-Boltzmann constant 
uj = Albedo of scattering 
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1 Introduction 


High porosity metallic foams are being studied as possible components of aero-space thermal protection 
systems. The performance of these systems needs to be characterized over a wide range of temperatures 
and pressures during ascent and re-entry. 1 In the future when such materials are specified as part of a 
vehicle program, part of the procurement process will involve certification that the material meets the 
specifications. 

To date there has been little research on experiment design for thermal characterization of high-porosity 
materials. The present research is intended to close this gap in the procurement cycle by investigating 
transient methods for accurate measurement of thermal properties in metal-foam materials. In this paper 
only numerical simulations are presented. 

A review of the pertinent literature is given next in the areas of optimal experiment design and heat 
transfer in high-porosity materials. Transient experiments combined with parameter estimation have been 
used for obtaining thermal properties for many years. 2 In these methods the desired parameters are found 
by non-linear regression between the experimental data (temperatures in this case) and a computational 
model of the experiment. Parameter estimation concepts have recently been applied to optimal experiment 
design for thermal characterization of uniform materials 3 and for materials with temperature- varying 
properties.' 1 One of us (Cole) has previously studied optimal experiment design for measurement of thermal 
conductivity in low-conductivity materials 5 and in layered materials. 6 The work was limited to a small 
rise in temperature so that radiation heat transfer was negligible. An optimality criterion was used to find 
the best experimental conditions for estimation of thermal properties. 

There have been many studies of glass-fiber insulation in which combined radiation and conduction 
heat transfer are present. Yuen et. al 9 used a detailed radiation model and measured optical properties to 
simulate the temperature response of a fibrous insulation material from first principles. The model involves 
a detailed determination of wavelength-dependent radiation exchange between multiple zones through the 
thickness of the material. Their results compared favorably to transient experimental tests carried out in 
1974 on LI9000 sintered silica insulation (rigid space shuttle tile). 

Insulation models based on optical properties, although science-based, do not provide direct evidence 
of the insulation’s ability to withstand reentry heating. Simplified insulation models, combined with direct 
tests, can provide compelling evidence that a thermal-protection system is ready for use on a human-piloted 
vehicle. Daryabeigi 7 used steady experiments to characterize low- density fibrous insulation over a wide 
temperature and pressure range and a transient experiment was used to simulate re-entry aerodynamic 
heating conditions. A combined radiation /conduction model was used, with the radiation described by 
a two-flux method assuming anisotropic scattering and a gray medium. Two models were discussed for 
combining solid and gas conduction into an effective conductivity, a parallel model and a model based 
on fiber orientation. Thermal properties estimated with both models yielded similar results, so it was 
not possible to choose one conduction model over the other. Low Rayleigh numbers and experimental 
verification ruled out free-convection heat transfer in the pores. 

There arc few studies of open-cell metallic foams. Zhu et. al 10 carried out simulations of a titanium 
foam material for the purpose of finding the minimum weight of the thermal protection by varying the pore 
size of the foam across its thickness. Their steady-state model included effective conductivity and radiation 
parameters. Compared to uniform materials, a density-graded material provides the same protection with 
less weight, with greater percentage savings for thinner insulation. 

Sullins and Daryabeigi 11 studied a Nickel-foam material over a wide range of temperatures and pressures 
with steady-state experiments combined with parameter estimation methods. The model for the heat 
transfer included a two-flux model for the radiation and combined gas/solid conduction. The effective 
conductivity was described by parallel gas-solid conduction, temperature-dependent gas conduction, and 
an efficiency factor for solid conductivity (as a small fraction of the bulk- metal conductivity). A gas-solid 
coupling term was also added to take into account an observed increase in conductivity at higher pressures 
and temperatures. The measurement of the model parameters was carried out with a non-linear regression 
technique, and the parameters reported. 

In this paper the methods of optimal experiment design are applied to a high-porosity nickel foam 
material, with thermal property values taken from Sullins and Daryabeigi. 11 Both conduction and radiation 
heat transfer are included. Based on a one-dimensional model of transient heat transfer, a large number 
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of simulated experiments have been studied to determine which experimental conditions provide the best 
estimates of the thermal properties. No laboratory experiments are reported. 

The paper is divided into several sections, as follows: analytical model; numerical solution; simulated 
experiments; optimality criterion; optimal experiment design; estimation of parameters; and conclusions. 


2 Analytical Model 


In this section the heat transfer model of Sullins and Daryabeigi 11 ’ 12 is used to model metallic foams. 
It is a one dimensional, two- flux model taking into account both conduction and radiation heat transfer. 
Contained within the model are five parameters considered to be intrinsic properties of the material that 
must be determined empirically. Three of these intrinsic properties are related to the radiative heat transfer 
and two are related to heat transfer by conduction. 

The transient heat equation with radiation is given by 


P dt dx dx ’ 


( 1 ) 


where T is the temperature, q" is the conduction heat flux, q" is the radiant heat flux, p is the density, 
and c is the specific heat. The conduction heat flux is given by 



(2) 


where k is the thermal conductivity. 

The gradient of the radiant heat flux is given by 

^=0{l- u )(G-4aT*), (3) 

where G is the incident radiation, w is the albedo of scattering, 0 is the extinction coefficient, and a is the 
Stefan-Boltzmann constant. The albedo of scattering is to be determined experimentally. The incident 
radiation is related to the radiant heat flux by 


,,__±dG 
Qr 3(3 dx ' 


( 4 ) 


Thus the heat equation (1) becomes 


pc 


dT 

dt 


dx \ dx ) 30 dx 2 


The incident radiation, G , satisfies the second order differential equation 

1 d 2 G 


30 2 (l - to) dx 2 

on the interior. On the boundary, G must satisfy: 

2 


+ G — 4<tT a 


( 5 ) 


(6) 


30(A) 

2 

30(a) 



= 4aT£, x = 0, 

(7) 

an 


= 4oT*, x = L, 

( 8 ) 


where is the temperature of the heater, T c is the temperature of the cold plate (300.6 K), — 0.85 is 

the emittance of the septum plate, and €2 = 0.92 is the emittance of the cold plate. 
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Next the material properties will be discussed for a specific high porosity nickel foam in a N 2 atmo- 
sphere. The product pc is the weighted average of the values for the metal (nickel) and the gas (nitrogen). 
Explicitly, 


pc = Pm Cm ( 1 - e) + PgCgt (9) 

The density, p g , and specific heat, c 5 , of the gas are temperature dependent and given in the appendix. 
For nickel, p m is taken to be 8900 and Cm = 444.0. The porosity, e, is 0.968 for nickel foam. 

Heat conduction is assumed to occur in parallel in the gas and solid with some coupling in the form 

k = ek g - h (1 - e) Fk^ + a ( Fk g k m ) 2 (10) 


where F is the efficiency and a is the coupling coefficient. The thermal conductivity of the metal is a 
function of temperature and is interpolated from tabulated data. 13 For the gas 


k 


9 ~ 


k* 

Z’ 


( 11 ) 


where k* is the temperature dependent conductivity of the gas at atmospheric pressure (see Appendix) 
and Z is given by 


Z = $ + 2$ 


2(2 — a)^Kn 
0(7 + l)Pr 


( 12 ) 


In this expression, the thermal accommodation coefficient, a, is taken to be one, 7 = 1.4 is the specific 
heat ratio, Pr is the temperature dependent Prandtl number (see Appendix), and Kn is the Knudson 
number. The Knudsen number is given by 


Kn = ^ (13) 

0 

where A is the gas molecular mean free path and S is the pore size of the metallic foam. The quantities 4> 
and were defined 11 as 


•©1 

11 

CD m 

Kn < 10 
Kn > 10 

(14) 

* - {5 

Kn < .01 
Kn > .01 ' 

(15) 


This was to provide an approximation of Z in the different ranges by neglecting one term when it was much 
smaller than the other. However, determining the optimum experiment requires finding the derivative of 
the temperature at the sensor with respect to the empirical parameters. Having a discontinuous expression 
for the gas conductivity resulted in inaccurate results for these derivatives when the simulated experiments 
were done near the boundary between two of the ranges for Kn. Therefore in this work 4> = = 1. 

The pore size for the metallic foam was assumed to be in the same form as given by Daryabeigi 12 and 
is given by 
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7T Df 

4 


(16) 


In this expression Df = \AE — 5 m is the diameter of the strut. The mean free path is given by 

KbT 
sj 2nPd 2 g 


(17) 


in which Kb = 1.38£ — 23 is the Boltzmann constant, P is the pressure, and d g = 3.798£ — 10 m is the 
gas collision diameter. 

The extinction coefficient is given by 


P = e- p 


(18) 
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where p is the density of the foam and e is the specific extinction coefficient. The specific extinction 
coefficient is taken to be a linear function of the temperature. Explicitly, 

e = eo 4- e\ T (19) 

The quantities eo, and ei are considered to be intrinsic properties of the media to be determined experi- 
mentally. 

Thus, for this model, heat transfer in a metallic foam is determined by the intrinsic properties F, a, 
eo, ei and uj. The goal is to find the optimal transient experiment to determine these parameters. The 
variables considered for different experiments are the heating power, sensor location, heating time and 
pressure. 


3 Numerical Solution 


The heat equation (5) and the radiation equation (6) are solved numerically. For each time step of the heat 
equation, the radiation equation (6) is discretized and solved. The solution is then used in a discretized 
form of the heat equation (5) to determine the temperatures for the next time step. The temperatures are 
then replaced with the new temperatures and the process is repeated. 

Discretizing (6) through (8) gives the tri-diagonal system of equations 
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where G{ and 7* are the radiation and temperature, respectively, at node i. The elements of the coefficient 
matrix are given by 


Ai 

Ci 

Bi 


2 

3/3? (1— w)Axi_ifAx7~i VAxi) ’ 
2(2 l-ca) 

3ftt2Asi_i ’ 


z = 2, . . . ,7V — 1, 
i = N, 


2(2— ex ) 

3/3eiA Xi » 

2 

3^^(T— LLf)Ax»(Az,:_i -{-Axi) 1 


i = 1, 

z = 2, . . . , TV — 1, 


1 - Ci, i = 1, 

1 -Ci-Ai, i = iV — 1, 
1 - A u i = N. 


( 21 ) 

( 22 ) 

(23) 


In discretized form the heat equation (5) is 


p c j _ i 


Ax, 


■ + pcy 


Ax 
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r pn -\- 1 rpri 

z± _i_ 

At 
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(24) 


3 Pj Ax j 


GJVO- 


where superscripts have been added to T and G to indicate the time step. Solving for TJ' +1 gives 


2 At 


'j-m-fl rpn 

J J AXj_i + pCjAx. 




(GJ-.-GJ)' 


( 25 ) 
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3/3 7 — i Axj _ i 
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The material properties of the metallic foam, pc, k , and 0, are all functions of temperature given by 
either polynomial approximations or interpolated from tabulated data. The temperature used for £* — piCi , 
ki , or 0i is the average temperature of the two endpoints of the interval Ax*. Explicitly for £ = £(T) 


Ci = C 


Ti + Ti 


i+i 


( 26 ) 


4 Simulated Experiments 


In this section the simulated experiment for determining the thermal properties of the metallic foam is 
described. The simulated experiment consists of heating the (one-dimensional) sample of nickel foam, at 
some pressure, for a particular time period with a heater of a specified power, while taking temperature 
measurements at one location in the body. Thus the variables of the experiment are the heater power, the 
sensor location, the pressure, the heating time, and the total duration of the experiment. 

Increasing the power to the heater increased the ability to determine the material properties. Since 
there would be physical limitations on the heater power supplied, rather than using the power of the heater 
as a variable between different experiments, an approximation of a heating apparatus used by Daryabeigi 14 
was made. This apparatus consisted of 2500 W cylindrical heating elements spaced 1.27 cm apart. This 
corresponds to about 310 MW/m 2 and an effective thickness of b — .00098 m. The heater and the sample 
were placed 5.08 cm apart. A cold plate, kept at room temperature, is considered to be along the other 
side of the sample. 

The temperature of the heater, node one in the discretized form of the heat equation, is determined 
from a lumped mass model. Under this model the heater must satisfy the differential equation 


PhChb 


dJ\ 

dt 


Qin Q.out i 


( 27 ) 


where p/ } , — 2600 kg/m 3 is the density of the heater and Ch = 800 J/(kg K) is the specific heat of the heater. 
The power supplied to the heater is Qi n — 310 MW/m 2 when the heater is powered and zero otherwise, while 
q oui is the heat loss to the sample and the surroundings. 

The losses from the heater include radiative transfer to the surroundings and to the sample. The 
radiative heat losses are given by 


<iout = e h a(Tjt - T 4 ) + €l a(T 4 - T 4 ), 


( 28 ) 


where th = .05 is the emittance of a water-cooled mirrored surface behind the heater, T 0 is the mirror 
temperature (300.6 K), and T s is the temperature of the surface of the sample. Gas conduction from the 
heater to the sample is also included, described by 

OT 

Qout = kg ■ ( 29 ) 

The total heat flow out of the heater is q out = q r out + q^uv 

A typical experiment starts at room temperature, then the heater is powered for a period of time while 
sampling the temperature at a particular location. Sampling continues for an additional period after the 
heater power is off. The heater model is incorporated into the transient finite difference code after the 
temperatures for each time step are found. 


5 Optimality Criterion 

The goal of this study is to examine a wide range of simulated experiments and determine the best 
experiment for determining the properties of the metallic foam. The suitability of individual experiments 
is determined from the sensitivity of the temperature with respect to the sought-after thermal properties 
of the sample, F, a, eo, ei and u j. Relabelling the parameters 6^, k = 1, . . . , 5, the normalized sensitivity 
coefficient for the kth. parameter, at the ith time step, is defined as 

X k (i) = b k ^±. ( 30 ) 
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Perturbing the parameters by a small amount, (1 + 8)b k , the sensitivity coefficients can be approximated 
by 


X k (i) a 


[Ti((l + 6)b k )-Ti(b k )] , 

- ,« 


1 S 


(31) 


The value of 5 — .0001 was found to give well-behaved values for X. 

In general an experiment is better if the sensitivity coefficients are larger. In addition, for experiments 
with more than one parameter the sensitivity coefficients must be linearly independent. For these reasons 
the optimality criterion used is 


D = 


|X T X| 

(nTl ax )p- 


(32) 


The sensitivity matrix, X, is defined by 


X T = 


' dr u dT , dT~ 
bl db 1 ,b2 db 1 ’"' ,b5 db 5 _ 


(33) 


where T is a vector of the temperatures at each sampling time. In (32), D is normalized by the maximum 
temperature, T max , squared and the number of sample temperatures, n , all taken to the power p (the 
number of parameters) in order to obtain a fair comparison among different experiments. 3 

The sensitivities are partial derivatives at a particular value of each parameter. The parameter values 
used for this study are given in Table 1, as reported by Sullins. 11 For all of the examples the thickness of 
nickel foam was 13.6 mm. 


Table 1: Parameter values for simulated experiments 


F a [(m K/W) 3 ] 

u e 0 [m 2 /kg] e 1 [m 2 /(kgK)] 

6.85E-3 3.89E2 

9.85 -2.63E-3 9.93E-1 


6 Optimal Experiment Design 

Our experience with thermal conductivity measurements, with no radiation present, suggested that it is 
best to heat the sample as rapidly as possible for a period of time while taking measurements at the 
surface closest to the heater. This is also the best approach for radiation/conduction materials. A variety 
of simulated experiments were explored by varying the pressure, the duration of the heating and the total 
duration of the experiment. Initially the simulated experiments were terminated at the time corresponding 
to the peak value of the optimality criterion, D. It was found, however, that extending the experiment 
duration increased the accuracy of the parameter estimation. For experiments at lower pressures radiative 
heat transfer is more dominant, while at higher pressures heat transfer due to conduction is more significant. 
It was thus expected that the optimum case would occur at moderate pressures since three of the parameters 
are related to radiation and two are related to conduction. 

The approximate conditions for the optimum experiment were found first by trial and error to find a 
rough value of the maximum of D, after which an optimization routine was run with the rough value as 
the initial guess. The maximum value of D found was 3.38E-16. This value resulted from an experiment 
conducted at 0.190 mm Hg in which the heater was run for 269 seconds and the temperature was sampled 
every 10 seconds for a total duration of 840 seconds. In Figure 1 the optimality criteria is graphed versus 
time for a few simulated experiments in which the optimum heating time was used and the pressure was 
varied around the optimal pressure. In Figure 2, D is graphed versus time with the pressure fixed at 0.19 
mm Hg for various heating times. As can be seen from the figure, the exact heating duration is not too 
critical. For heating times between 260 and 290 seconds the maximum value of D only varies by about four 
percent. In general if the sample is continuously heated the D value rises to a maximum and then begin to 
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decrease as the system reaches steady state. If instead the heater is shut off, at some time, the the value 
of D dramatically increases until the system again nears steady state (room temperature). The rise in the 
D values when the heater is shut off is so dramatic that in most cases it is impossible to see the rise in the 
D values during heating if the entire experiment is graphed on the same scale. In Figure 3 the same curves 
given in Figure 2 are repeated along with a curve corresponding to constant heating for 1200 seconds. 
The separation in the four curves prior to the heater being shut off is due to each experiment reaching a 
different maximum temperature which is included in the normalization. The best time for heating (269 
seconds) appears to be near the location where the D curve changes from concave up to concave down. In 
Table 2 the maximum value of D obtained for various pressures and heating times is summarized. 


Table 2: Maximum D values for experiments 


P [mmHg] 



Heating Duration [sec] 




250 

260 

269 

280 

290 

300 

0.09 

1.30E-16 

1.30E-16 

1.30E-16 

1.30E-16 

1.28E-16 

1.27E-16 

0.14 

2.69E-16 

2.71E-16 

2.72E-16 

2.71E-16 

2.69E-16 

2.67E-16 

0.19 

3.19E-16 

3.23E-16 

3.38E-16 

3.25E-16 

3.24E-16 

3.22E-16 

0.24 

4.32E-16 

4.11E-16 

3.93E-16 

3.78E-16 

3.65E-16 

3.53E-16 

0.29 

4.73E-16 

4.51E-16 

4.32E-16 

4.15E-16 

4.01E-16 

3.89E-16 


It should be noted that for the optimal experiment one would want to stop taking data as soon as D 
reached its maximum value. From Figure 1 or Figure 2 this is at 840 seconds. This time is also not too 
critical, since for any value between 830 seconds and 900 seconds D will still be within one percent of its 
maximum value. Thus for this model the optimal experiment for determining the material properties is 
for the sample to be heated for 270 seconds at a pressure of 0.19 mm Hg with sampling every 10 seconds 
until 840 seconds have elapsed. 

In Figure 4 the heat flux for the optimal experiment is given as a function of time. The total heat flux 
is approximately a step function while the radiative heat flux and heat flux due to conduction are of the 
same order, as expected. The temperature profile for the optimal experiment is given in Figure 5. One 
can see from the temperature profile the system has not reached steady state when the heater is shut off. 


7 Estimation of Parameters 

In this section the parameters are determined using three simulated experiments. Experiment 1 is the 
experiment found using the optimality criterion in the previous section. Experiment 2 is the same as 
experiment 1 except data is taken for 2000 seconds. Experiment 3 consists of two heating events concate- 
nated into one simulated experiment, where one of the heating events is dominated by conduction flux 
and the other is dominated by radiative flux. Data for all of the simulated experiments are obtained by 
calculating the exact temperatures using the values of the parameters given in Table 1 and then corrupt- 
ing them with random noise with a standard deviation of 0.08 K. The standard deviation was chosen to 
correspond with the maximum standard deviations reported by Sullins 11 for actual experiments in the 
same temperature ranges. A nonlinear Levnberg-Marquardt parameter estimation routine is used to fit 
the model parameters. 15 

For simulated experiment 1 the optimum experiment determined in the previous section was used. It 
consists of running the heater for 269 seconds and taking measurements until 840 seconds elapsed. The 
gas pressure is 0.19 mm Hg. The residuals between the exact and perturbed (noise added) temperatures, 
Tpert — Texact , and between the exact temperatures and those calculated using the estimated parameters, 
Test—Texact , are given in Figure 6. The residuals between the exact temperatures and those calculated with 
the estimated parameters is less than 0.5 K and is biased towards overestimating the exact temperatures. 
In addition the error in the calculated temperatures exceeds the random noise added to the exact values. 
The estimated values of the parameters along with their 90% confidence intervals are given in Table 3. 

Since at 840 seconds the temperature was still significantly above room temperature, for simulated 
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Table 3: Parameter values for simulated experiment 1 



eo 

ei 

D 

F 

a 

Exact 

9.850E+00 

-2.630E-03 

9.930E-01 

6.850E-03 

3.890E+02 

Estimated 

9.838E+00 

-2.619E-03 

9.930E-01 

6.582E-03 

4.696E-F02 

Percent Error 

-0.127 

-4.17E-01 

4.38E-03 

-3.91E00 

2.07E01 

90 % Confidence Interval 

6.551E-03 

4.922E-06 

3.152E-05 

7.956E-05 

2.821E+01 


experiment 2 the same heating time was used, but the temperature sampling was extended to 2000 seconds. 
By 2000 seconds the temperature at the sensor had cooled to 305 K. The results were much better both 
in terms of the residuals (See Figure 7) and the parameter estimation (see Table 4). The errors in the 
temperatures were less than .05 K while the parameters were found accurately within about one percent. 


Table 4: Parameter values for simulated experiment 2 



eo 

ei 

UJ 

F 

a 

Exact 

9.850E+00 

-2.630E-03 

9.930E-01 

6.850E-03 

3.890E-F02 

Estimated 

9.853E+00 

-2.633E-03 

9.930E-01 

6.874E-03 

3.819E+02 

Percent Error 

3.489E-02 

1.225E-01 

-1.233E-04 

3.521E-01 

-1.824E+00 

90 % Confidence Interval 

6.569E-03 

4.925E-06 

3.191E-05 

5.118E-05 

1.529E+01 


Table 4 shows that simulated experiment 2 estimated all of the parameters accurately, except for the 
coupling coefficient, a. Various other experiments with single heating events were tried, without success, in 
order to find the coupling coefficient more accurately. The coupling coefficient is only significant when the 
conduction through the gas is significant. This means a will be more important in experiments conducted 
at high pressure. In addition it should be easier to determine a from an experiment in which the radiative 
heat transfer is less important. Based on this information, simulated experiment 3 was constructed from 
a pair of heating events: one that was optimal for determining the conduction parameters; and one that 
was optimal for finding the radiation parameters. Noise was added to the data from this pair of heating 
events. The data sets were then concatenated as one data set for the parameter estimation. It should be 
noted that although optimal heating events were found for both the conduction and radiation parameters 
with the other parameters fixed, for the analysis of the concatenated data, all of the parameters were 
considered to be unknown. 

For the conduction-optimal heating event, increasing pressure was found to improve the estimation 
of a and F. However, an increase in pressure from 1 to 10 atmospheres increased the D values by less 
than one percent. Since it would be much easier to obtain a pressure of one atmosphere this was used 
for the conduction-optimum heating event. In this case, it was also found that using more power in the 
heater was better. Therefore the only variables to determine was how long to power the heater and how 
long to sample the temperature. From Figure 4 it can be seen that the flux due to conduction reaches 
a peak at about 20 seconds. At one atmosphere this peak moves to about 23.6 seconds, which was used 
for the conduction-optimal heating event. The temperature was sampled every second instead of every 
ten seconds. It was found that a data duration of 150 seconds provided good estimation of conduction 
parameters. Strict use of the D parameter would have resulted in a shorter experiment, but as with the 
single heating event the sample had not cooled very close to steady state by that time. 

For the radiation-optimal heating event, D is the largest for a vacuum (no gas conduction), Sullins 11 
reported using a minimum pressure of 0.0001 mm Hg so this value was chosen as a practically obtainable 
minimum pressure. For the experiments discussed so far, the initial heating had a significant conduction 
flux (see Figure 4). In Figure 8, the ratio of the conduction heat flux to the radiation heat flux is plotted 
for a few different heating times. It can be seen from the figure that if the sample is heated for 30 seconds 
this ratio is nearly zero for a large portion of the experiment. This corresponds to the experiment being 
dominated by radiation. For this reason the radiation-optimal experiment was to heat the sample for 30 
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seconds and then data was sampled for 100 seconds. Again a longer sampling time was used than the D 
value indicated was optimum. 

Simulated experiment 3 consisted of the conduction-optimal and radiation-optimal heating events com- 
bined together. The residuals for simulated experiment 3 are given in Figure 9. The estimated parameters 
and their 90% confidence intervals are given in Table 5. As can be seen in the figure and the table, the 
temperatures and the parameters are estimated accurately. 


Table 5: Parameter values for combined conduction and radiation experiments 



eo 

ei 

u> 

F 

a 

Exact 

9.850E+00 

-2.630E-03 

9.930E-01 

6.850E-03 

3.890E+02 

Estimated 

9.844E+00 

2.626E-03 

9.930E-01 

6.838E-03 

3.905E+02 

Percent Error 

6.271E-02 

1.666E-01 

-2.696E-05 

1.792E-01 

-3.960E-01 

90 % Confidence Interval 

1.944E-02 

1.668E-05 

2.818E-05 

9.107E-06 

8.771E-01 


The sensitivity coefficients (normalized by their maximum value) are plotted for simulated experiment 
1 in Figure 10. The sensitivity for a has nearly the same shape as the sensitivity for F. This indicates the 
sensitivities are almost linearly independent and is probably the reason it is hard to determine the value 
for a. A different choice of conduction parameters could probably be found that could be determined more 
accurately from simulated or actual experiments. 

Simulated experiment 3 provides more accurate results than either experiment 1 or experiment 2. For 
all the experiments, the strict use of the normalization of D by the number of samples taken tended to end 
the experiments before all useful data was taken. Unnormalized values of X T X values would indicate the 
sample should be heated to steady-state and then cooled back to room temperature, which would require a 
very long experiment. The best data duration appears to be somewhere between the duration determined 
by strict use of the maximum value of the normalized and unnormalized values of X T X. 

8 Conclusion 

For future use of high porosity materials in thermal protection systems for reusable launch vehicles, their 
performance needs to be characterized over a wide range of temperatures and pressures. The goal of 
this research is to develop general-purpose protocols for designing transient experiments to measure ther- 
mal properties of high porosity materials. As a first step, a metal foam material was studied with an 
experimentally validated thermal model. 

The D optimality criterion and consideration of the heat fluxes can lead to suitable experiments for 
determining the properties. Maximization of the optimality criterion D lead to an experiment containing 
one heating event that could be used to determine the radiation and conduction parameters. Strict use 
of the D parameter leads to experiments that were too short in the sense that more accurate parameter 
estimates could be obtained from simulated data sets of longer duration. By considering the model itself 
and in particular its conduction and radiation behavior, a set of two shorter heating events were found 
that provided more accurate parameter estimates. 

The radiation parameter u was easily found to a high degree of accuracy by all of the experiments. 
The coupling coefficient, a, was the hardest parameter to determine. The sensitivity of the of the coupling 
coefficient and the sensitivity of the solid conduction efficiency, F, are similar in shape which indicates 
they are close to linearly dependent. The coupling coefficient was an ad-hoc parameter added to Sulims’ 
model in order to obtain agreement with experimental data at high pressures. This suggests that there 
may be a different pair of conduction parameters that could adequately model the conduction process and 
provide more robust estimation than exhibited here. 
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Appendix 

The temperature dependent properties for nitrogen are from Daryabeigi. 12 

8.1 Temperature Dependent Nitrogen Properties 


p g = -1.734 E - 4 + 342.216/T (34) 

Cg = 1.083.545 - 0.328T + 6.949E - 4T 2 - 2.82 E - 7 T 3 (35) 

k g = 2.048E - 3 + 8.751F - 5T - 2.462E - 8T 2 (36) 

Pr = 0.854 - 7.085F - 4T 4- 9.008F - 7 T 2 - 3.207P - 10T 3 (37) 


References 

1 Dorsey, J. T. et al., “Metallic Thermal Protection System Technology Development: Concepts, Re- 
quirements and Assessment Overview”, Proceedings , 40th Aerospace Sciences Meeting , Jan 14-17, 2002, 
Reno, Nevada, paper AIAA 2002-0502. 

2 Beck, J. V. and Arnold, K. J., Parameter Estimation , 1977, Wiley, NY. 

3 Taktak, R., Beck, J. V., and Scott, E. P., “Optimal experimental design for estimating thermal prop- 
erties of composite materials,” International Journal of Heat and Mass Transfer , Vol. 36, No. 12, 1993, 
pp. 2977-2986. 

4 Dowding, K. J. and Blackwell, B. F., “Sensitivity Amalysis for Nonlinear Heat Conduction,” Journal 
Heat Transfer , 123, 1-10 (2001). 

5 Cole, K. D., “Thermal Characterization of Functionally Graded Materials-Design of Optimal Exper- 
iments,” 8th Joint AIAA/ ASME Thermophysics and Heat Transfer Conference , St. Louis, MO, June 
24-27, 2002, paper AIAA 20002-2882. 

6 Cole, K. D., “Analysis of Photothermal Characterization of Layered Materials-Design of Optimal Ex- 
periments,” 15th Symposium on Thermophysical Properties, Boulder, CO, June 22-27, 2003. 

7 Daryabeigi, K., “Heat Transfer in High- Temperature Fibrous Insulation,” AIAA J . Thermophysics and 
Heat Transfer, vol. 17, no. 1, 2003, pp. 10-20. 

8 Bhattacharyya, R., “Heat Transfer Model for Fibrous Insulations,” Thermal Insulation Performance , 
D. L. McElroy and R. P Tye, editors, American Society for Testing and Materials Technical Publication 
718, Philadelphia, pp. 272-286, (1980). 

9 Yuen, W. W., Takara, E., and Cunnington, G., “Combined Conductive/Radiative Heat Transfer in 
High Porosity Fibrous Insulation Materials: Theory and Experiment,” Proceedings , 6th ASME-JSME 
Joint Thermal Engineering Conference , March 16-20, 2003, paper TED-AJ03-126. 

10 Zhu H., et. al “Minumum Mass Design of Insulation Made from Functionally Graded Design”, Proceed- 
ings , 43rd AIAA Structures, Structural Dynamics, and Materials Conference, Denver CO, April 22-25, 
2002, paper AIAA-2002-1425. 

11 Sullins, A. D., and Daryabeigi, K., “Effective Thermal Conductivity of High Porosity Open Cell Nickel 
Foam, ” Proceeding of the 35th AIAA Thermophysics Conference, Anaheim CA, June 11-14, 2001, 
paper AIAA 2001-2819. 

12 Daryabeigi, K., Design of High Temperature Multilayer Insulation for Reusable Launch Vehicles, Ph.D. 
dissertation, University of Virginia, May 2000. 


12 



13 Rohsenow, W. M., Hartnett, J.P., and Cho, Y.I., Handbook of Heat Transfer, 3 rd Edition , The McGraw- 
Hill Companies, Inc., 1998. 

14 Daryabeigi, K., Knutson, J. R., and Sikora, J. G., “Thermal Vacuum Facility for Testing Thermal 
Protection Systems,” NASA TM-2002-211734, June 2002. 

15 Press, W. H., Numerical Recipes in FORTRAN 77, Cambridge University Press, New York, 1996. 



Figure 1: D for 269 seconds of heating at various pressures in mm Hg 
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Figure 2: D for P=0.19 mm Hg and various heating times 


1000 


1200 



14 



temperature (K) 







Residu 





Residuals (K) 



Figure 8: Ratio of conduction to radiation flux for different heating times in seconds. 



Figure 9: Temperature residuals for simulated experiment 3. 
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Figure 10: Sensitivities for conduction portion of simulated experiment 3 
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